/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "linal/common.hxx"
#include "linal/util.hxx"

#include "linal/format.hxx"

#include "linal/operation.hxx"

#include "linal/device/host/operation.hxx"

namespace linal {

  static int apply_q_inc_lhfc(Hier_ A, Hier_ T, Hier_ W_1, Hier_ B);
  static int apply_q_inc_lnfc(Hier_ A, Hier_ T, Hier_ W_1, Hier_ B);

  /*!
    left, conj_transpose, forward, columnwise
  */
  int apply_q_inc_lhfc(Hier_ A, Hier_ T, Hier_ W_1, Hier_ B) {
      
    FLA_Obj ATL,   ATR,      A00, A01, A02, 
            ABL,   ABR,      A10, A11, A12,
                             A20, A21, A22;

    FLA_Obj TTL,   WTR,      T00, W01, W02, 
            TBL,   TBR,      T10, T11, W12,
                            T20, T21, T22;

    FLA_Obj BT,              B0,
            BB,              B1,
                             B2;

    FLA_Part_2x2( ~A, &ATL, &ATR,
                      &ABL, &ABR, 0, 0, FLA_TL );
    FLA_Part_2x2( ~T, &TTL, &WTR,
                      &TBL, &TBR, 0, 0, FLA_TL );
    FLA_Part_2x1( ~B, &BT,
                      &BB,        0, FLA_TOP );

    while ( FLA_Obj_min_dim( ABR ) > 0 ){

      FLA_Repart_2x2_to_3x3( ATL,  ATR,       &A00,  &A01, &A02,
                                              &A10,  &A11, &A12,
                             ABL,  ABR,       &A20,  &A21, &A22,
                             1, 1, FLA_BR );
      FLA_Repart_2x2_to_3x3( TTL,  WTR,       &T00,  &W01, &W02,
                                              &T10,  &T11, &W12,
                             TBL,  TBR,       &T20,  &T21, &T22,
                             1, 1, FLA_BR );
      FLA_Repart_2x1_to_3x1( BT,                &B0,
                                                &B1,
                             BB,                &B2,        1, FLA_BOTTOM );
      //------------------------------------------------------------
      Hier_ 
        A_11(A11), A_21(A21), 
        T_11(T11), T_21(T21),
        B_1(B1), B_2(B2);

      void
        *buf_A_21 = A_21.get_buffer(),  *buf_T_21 = T_21.get_buffer(),
        *buf_W_1  = W_1.get_buffer(),   *buf_B_1  = B_1.get_buffer(),
        *buf_B_2  = B_2.get_buffer();
      

      int
        cs_A_21 = A_21.get_cs(),  cs_T_21 = T_21.get_cs(),
        cs_W_1  = W_1.get_cs(),   cs_B_1  = B_1.get_cs(),
        cs_B_2  = B_2.get_cs();

      apply_q( FLA_LEFT, FLA_CONJ_TRANSPOSE, 
               FLA_FORWARD, FLA_COLUMNWISE,
               A_11, T_11, W_1, B_1 );

      // do not change the order of loop
      for (int i=0;i<B_2.get_m();++i) {                 
        for (int j=0;j<B_2.get_n();++j) {

#pragma omp task firstprivate(A_21, T_21, W_1, B_1, B_2, i, j) 
          FLA_Apply_Q2_UT_internal( FLA_LEFT, FLA_CONJ_TRANSPOSE, 
                                    FLA_FORWARD, FLA_COLUMNWISE,
                                    LINAL_ELEM_FLA(_A_21, i,0), 
                                    LINAL_ELEM_FLA(_T_21, i,0),
                                    LINAL_ELEM_FLA(_W_1, 0,j),  
                                    LINAL_ELEM_FLA(_B_1, 0,j),
                                    LINAL_ELEM_FLA(_B_2, i,j),
                                    fla_apq2ut_cntl_leaf );

        }

#pragma omp taskwait	  

      }

      //------------------------------------------------------------

      FLA_Cont_with_3x3_to_2x2( &ATL,  &ATR,       A00, A01,  A02,
                                                   A10, A11,  A12,
                                &ABL,  &ABR,       A20, A21,  A22,
                                FLA_TL );
      FLA_Cont_with_3x3_to_2x2( &TTL,  &WTR,       T00, W01,  W02,
                                                   T10, T11,  W12,
                                &TBL,  &TBR,       T20, T21,  T22,
                                FLA_TL );
      FLA_Cont_with_3x1_to_2x1( &BT,                B0,
                                                    B1,
                                &BB,                B2,     FLA_TOP );
    }
    return true;
  }

  /*!
   left, no_transpose, forward, columnwise
  */
  int apply_q_inc_lnfc(Hier_ A, Hier_ T, Hier_ W_1, Hier_ B) {

    FLA_Obj ATL,   ATR,      A00, A01, A02, 
            ABL,   ABR,      A10, A11, A12,
                             A20, A21, A22;

    FLA_Obj TTL,   WTR,      T00, W01, W02, 
            TBL,   TBR,      T10, T11, W12,
                             T20, T21, T22;

    FLA_Obj BT,              B0,
            BB,              B1,
                             B2;

    FLA_Part_2x2( ~A, &ATL, &ATR,
                      &ABL, &ABR, 0, 0, FLA_BR );
    FLA_Part_2x2( ~T, &TTL, &WTR,
                      &TBL, &TBR, 0, 0, FLA_BR );
    FLA_Part_2x1( ~B, &BT,
                      &BB,        0, FLA_BOTTOM );

    while ( FLA_Obj_min_dim( ATL ) > 0 ){

      FLA_Repart_2x2_to_3x3( ATL,  ATR,       &A00,  &A01, &A02,
                                              &A10,  &A11, &A12,
                             ABL,  ABR,       &A20,  &A21, &A22,
                             1, 1, FLA_TL );
      FLA_Repart_2x2_to_3x3( TTL,  WTR,       &T00,  &W01, &W02,
                                              &T10,  &T11, &W12,
                             TBL,  TBR,       &T20,  &T21, &T22,
                             1, 1, FLA_TL );
      FLA_Repart_2x1_to_3x1( BT,              &B0,
                                              &B1,
                             BB,              &B2,        
                             1, FLA_TOP );
      //------------------------------------------------------------
      Hier_ 
        A_11(A11), A_21(A21), 
        T_11(T11), T_21(T21),
        B_1(B1), B_2(B2);

      void
        *buf_A_21 = A_21.get_buffer(),  *buf_T_21 = T_21.get_buffer(),
        *buf_W_1  = W_1.get_buffer(),   *buf_B_1  = B_1.get_buffer(),
        *buf_B_2  = B_2.get_buffer();


      int
        cs_A_21 = A_21.get_cs(),  cs_T_21 = T_21.get_cs(),
        cs_W_1  = W_1.get_cs(),   cs_B_1  = B_1.get_cs(),
        cs_B_2  = B_2.get_cs();

      // do not change the order of loop
      for (int i=(B_2.get_m() - 1);i>=0;--i) {
        for (int j=(B_2.get_n() - 1);j>=0;--j) {
            
#pragma omp task firstprivate(A_21, T_21, W_1, B_1, B_2, i, j) 
          FLA_Apply_Q2_UT_internal( FLA_LEFT, FLA_NO_TRANSPOSE, 
                                    FLA_FORWARD, FLA_COLUMNWISE,
                                    LINAL_ELEM_FLA(_A_21, i,0), 
                                    LINAL_ELEM_FLA(_T_21, i,0),
                                    LINAL_ELEM_FLA(_W_1, 0,j),
                                    LINAL_ELEM_FLA(_B_1, 0,j),
                                    LINAL_ELEM_FLA(_B_2, i,j),
                                    fla_apq2ut_cntl_leaf );

        }

#pragma omp taskwait	  

      }

      apply_q( FLA_LEFT, FLA_NO_TRANSPOSE, 
               FLA_FORWARD, FLA_COLUMNWISE,
               A_11, T_11, W_1, B_1 );

      //------------------------------------------------------------

      FLA_Cont_with_3x3_to_2x2( &ATL,  &ATR,       A00, A01,  A02,
                                                   A10, A11,  A12,
                                &ABL,  &ABR,       A20, A21,  A22,
                                FLA_BR );
      FLA_Cont_with_3x3_to_2x2( &TTL,  &WTR,       T00, W01,  W02,
                                                   T10, T11,  W12,
                                &TBL,  &TBR,       T20, T21,  T22,
                                FLA_BR );
      FLA_Cont_with_3x1_to_2x1( &BT,                B0,
                                                    B1,
                                &BB,                B2,     
                                FLA_BOTTOM );
    }
    return true;
  }


  int apply_q_inc(int side, int trans, int direct, int storev,
                  Hier_ A, Hier_ T, Hier_ W_1, Hier_ B) {

    if ( !B.get_m() || !B.get_n() ) return true;

    /*
    LINAL_ERROR_DETAIL( (A.get_m() == T.get_m()) &&
                        (A.get_n() == T.get_n()) ,
                        "%s\n dim A = ( %4d, %4d )\n dim T = ( %4d, %4d )\n",
                        LINAL_ERROR_MATRIX_DIMENSION_MISMATCH,
                        A.get_m(), A.get_n(),
                        T.get_m(), T.get_n() );
    */

    LINAL_ERROR_DETAIL( (A.get_m() == B.get_m()),
                        "%s\n length A = %4d\n length B = %4d\n",
                        LINAL_ERROR_MATRIX_DIMENSION_MISMATCH,
                        A.get_m(), B.get_m() );

    LINAL_ERROR_DETAIL( (W_1.get_n() == B.get_n()),
                        "%s\n width W_1 = %4d\n width B   = %4d\n",
                        LINAL_ERROR_MATRIX_DIMENSION_MISMATCH,
                        W_1.get_n(), B.get_n() );


    switch (side) {
    case FLA_LEFT:
      {
        switch (trans) {
        case FLA_NO_TRANSPOSE: 
        case FLA_CONJ_NO_TRANSPOSE: 
          apply_q_inc_lnfc( A, T, W_1, B ); 
          break;
        case FLA_TRANSPOSE:
        case FLA_CONJ_TRANSPOSE:    
          apply_q_inc_lhfc( A, T, W_1, B ); 
          break;
        }
        break;
      }
    case FLA_RIGHT:
      {
        LINAL_ERROR( false, LINAL_ERROR_NOT_YET_IMPLEMENTED );
        break;
      }
    }
    return true;
  }
}
